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Abstract 

This  project  concerns  the  development  of  an  image  reconstruction  algorithm 
which  removes  noise  from  a  corrupted  image  and  recovers  patterns  and  textures.  The 
textures  this  algorithm  detects  are  edges  and  smooth,  gray  scale  regions.  In  previous 
attempts  of  this  kind,  images  were  recovered  with  nice  edges,  but  blockiness  remained  in 
the  smooth  regions.  Our  goal  is  to  recover  both  types  of  image  composition.  We  do  this 
by  a  least  squares  minimization  method  and  a  regularization  process.  The  regularization 
is  the  key  element  that  allows  us  to  differentiate  the  edge  and  the  smooth  regions.  The 
computation  generated  by  the  algorithm  tends  to  create  a  large  system  to  solve.  We  will 
use  an  algebraic  multigrid  solver  to  minimize  the  time  for  solving  these  large  systems. 
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Introduction 


Signal  processing  and  specifically  image  restoration  span  a  wide  variety  of  uses. 
Whether  the  image  is  an  old  photograph,  an  infrared  scan,  or  a  satellite  image,  it  is  clear 
that  noise  filtering  and  shape  recovery  are  two  important  aspects  of  creating  a  clear 
picture. 

The  first  question  we  must  answer  is  what  is  the  need  for  this  type  of  image 
recovery.  Take  a  satellite  image  for  example.  Due  to  many  different  reasons,  a  satellite 
image  may  be  compromised  by  noise.  Whether  this  noise  is  physical  as  in  partial  cloud 
cover,  patterned  as  in  noise  added  during  the  transmission  of  the  signal,  or  random  due  to 
defects  in  equipment,  we  must  try  to  overcome  these  different  sources  of  noise.  But 
decreasing  noise  is  not  the  only  task.  Preserving  the  textures,  called  shape  preservation, 
is  equally  important.  A  noiseless  picture  is  not  much  benefit  unless  the  image  you  are 
observing  is  an  accurate  representation  of  what  is  being  observed. 

In  developing  an  algorithm,  it  is  beneficial  to  look  at  the  one-dimensional  case. 
Consider  a  one-dimensional  image,  I(x),  as  in  Fig  1.  Notice  the  three  different  types  of 


texture  regions  that  compose  this  picture.  There  are  three  flat  regions:  two  with  a  level  of 
zero,  and  one  with  a  level  of  1.  The  vertical  line  represents  an  edge,  while  the  sloped  line 
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segment  represents  a  smooth  region.  A  typical  noisy  observation  of  I(x)  might  occur  as 
what  is  seen  in  Fig  2.  Depending  on  the  amount  of  noise  in  the  image,  the  human  eye 


itself  cannot  perceive  the  desired  information.  That  is  why  we  need  a  method  to  recover 
the  image.  A  typical  reconstructed  image  is  shown  in  Fig  3. 


Notice  that  the  difference  between  a  reconstructed  edge  and  a  smooth  region  is  minor.  A 
good  algorithm  will  detect  the  difference  and  display  them  accordingly.  This  so-called 
shape  preserving  is  the  other  aspect  of  image  recovery  that  is  important. 

As  you  can  see  there  are  two  different  aspects  to  the  noisy  image.  One  is  the  high 
frequency  noise  that  the  noise  created.  The  other  is  a  low  frequency  effect  that  represents 
the  actual  shape  of  the  image.  First  we  describe  the  reconstruction  though  the  use  of 
partial  differential  equations.  We  apply  the  heat  equation,  which  is  well  known  to  have 
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smoothing  properties:  solving  the  heat  equation  with  a  discontinuous  initial  condition  will 
yield  a 


du  _  2 

dt  dx2 
u(0,x )  =  y(x) 


[1] 


continuous  and  smooth  solution.  This  smoothing  however,  tends  to  smear  the  edges. 

In  order  to  retain  and  recover  the  edges,  it  has  been  proposed  that  replacing  the 
heat  equation  with  the  non-linear  diffusion  equation. 


du 

¥ 


(  du'] 


G  \Ur 


d_ 

dx 


[2] 


V  J 

The  two  dimensional  version  of  this  nonlinear  diffusion  equation  is  known  as  the  mean 
curvature  flow.  The  nonlinear  diffusion  thus  provides  smoothing  in  the  direction  of  the 
tangent,  rather  than  in  the  normal  direction.  Thus,  the  edges  are  not  smeared  by  this 
procedure. 

Now  we  describe  an  alternative  approach  based  on  a  variational  formulation.  We 
consider  the  following  minimization 


J{u )  =  —  ||m(x)  -  y{x)^ dx  +  —  ^\ux^ dx  . 


[3] 


The  first  term  of  the  cost  functional,  J,  is  the  least  squares  fit,  u(x),  to  the  observed  image, 
y(x).  Minimizing  the  norm  of  the  difference  between  the  solution  and  the  observed 
image  insures  that  the  solution  stays  relatively  close  to  the  observed  image.  This  is  an 
important  term,  however  if  it  were  the  only  term,  the  minimum  solution  would  always  be 
the  observed  data.  The  second  term  is  called  Gaussian  Regularization.  This 
regularization  term  is  to  control  the  total  variation  of  the  reconstructed  image,  so  the 


4 


oscillation  of  the  constructed  image  is  controlled.  The  necessary  and  sufficient 
optimality  condition  is  given  by 

-  fiu„  +  u-y  =  0 .  [4] 

This  is  the  steady  state  of  the  forced,  time-dependent  heat  equation 

In  order  to  reconstruct  the  image  with  edges  we  consider  a  regularization  with  a  bounded 
variation  of  the  image,  u.  We  consider  instead  of  Eq.  3,  the  cost  functional 

J(u)=^\u(x)-y(x)^dx+ fi^\ux\dx  [6] 

It  will  be  shown  that  the  optimality  condition  of  this  problem  is  written  as 

lu  \ 

-Pri  +u  =  y  [7] 

{\Ux\  )x 

Then  we  can  observe  that  the  nonlinear  diffusion  term  is  similar  to  the  one  that  appears  in 
the  method  mentioned  above.  Now  it  is  well  known  that  this  BV  regularization  method 
works  well  for  images  with  edges,  but  our  objective  is  to  introduce  a  new  class  of 
regularization  to  avoid  the  blockiness  found  in  the  BV  regularization  method. 
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Regularization  Method 

In  this  section,  we  will  describe  our  idea  and  approach  using  the  one-dimensional 
case.  This  approach  can  then  be  extended  to  multiple  dimensions.  The  regularization 
method  we  employed  is  based  on  bounded  variations  of  the  image.  Here  our  domain,  Q,, 
is  the  interval,  (0,1).  The  function  we  would  like  to  minimize  is 

J(u)  =  —  J|w(jc)  -  zixtfdx  +  PxBV(u)  +  fi2BV(ux) ,  [8] 

^  o 

where  u  e  Hl0  (0,1) .  BV (u)  denotes  the  bounded  variation  semi-norm  of  u.  If  u  is 
absolutely  continuous,  then 

i 

BV(u)  =  j\ux\dx.  [9] 

o 

Thus  in  the  rest  of  our  discussion,  we  will  continue  to  use  the  right  hand  side  of  the 
previous  equation  as  our  notation  for  BV(u).  We  use  a  bounded  variational  approach 
since  this  allows  for  jump  discontinuities,  which  will  be  important  later  on  in  our 
discussion.  This  first  regularization  term,  BV(u),  represents  a  total  variation  of  the  image 
and  helps  to  restore  regions  with  edges.  In  other  words,  in  a  one-dimensional  image,  an 
image  that  is  a  piece-wise  constant  would  be  readily  restored. 

Now  while  this  term  alone  would  work  well  for  an  image  with  such  block 
characteristics,  it  does  poorly  in  restoring  what  would  be  represented  in  a  one¬ 
dimensional  image  as  a  slope.  Smooth  slopes  will  be  restored  to  a  series  of  jagged  stair 
steps.  Basically  it  will  try  to  give  edge-like  characteristics  to  the  smooth  region.  In  a 
gray  scale  image,  the  gray  would  be  restored  to  blocky  regions  of  black,  white,  and  gray 
levels  if  we  employ  only  the  first  BV  term.  To  restore  the  intermediate  gray  scale,  we 
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introduce  the  second  BV  term.  This  is  a  bounded  variation  on  the  derivative,  ux.  If  ux  is 
absolutely  continuous,  then 


BV(ux) 


1 


[10] 


where 

[w* ]($,.)  =  ux  (s* )  -  ux(s~ )  is  the  jump  discontinuity  of  the  derivative  at  x  =  s;.  This  BV 

term  is  the  shape  preserving  entity,  which  we  need  to  prevent  blockiness. 

Now  these  two  terms  can  be  used  to  relate  the  effects  on  each  other  on  the  way 
the  image  is  restored.  The  importance  of  the  relative  weight  factors,  Pi  and  p2,  can  be 
seen  by  letting  p2  be  zero  thereby  taking  away  gray  scale  resolution.  On  the  other  hand 
letting  Pi  equal  to  zero  would  make  edges  tend  to  disappear.  The  goal  is  to  find  a  balance 
between  these  two  extremes  to  find  a  way  to  reproduce  both  types  of  image  composition. 
For  example  we  can  take  a  piecewise  linear  function  on  the  interval  (0,1)  that  has  finite 
restoration  energy 

m  m-1 

Q(u)  =  frBV(u)  +  J32BV(u)  =  A  -*,•)  +  AEKu  “  4| .  [1 1] 

1=1  1=1 

where  u  is  linear  on  each  interval  (si,Si+i)  with  slope  X\. 

In  the  case  of  a  discontinuous  image  we  let  <j)  be  a  piecewise  constant  function, 


where  =  0  on  (0,  — )  and  <|)  =  1  on  ( — , 
2  2 


1)  can  be  approximated  by  the  piecewise  linear, 


continuous  function 
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<t>s  = 


IH 


+  - 


XE 


XE 


XE 


(o}-^ 

2  , 

1-5  1  +  6} 


fl+l^ 


[12] 


This  allows  us  to  bound  the  variation  of  ux,  for  otherwise  we  would  have  BV(ux)  = 


Now  we  have 


Q(<t>s)  =  fr  + 


2A 


[13] 


Ideally  we  would  like  the  restoration  energy  to  have  equal  contributions  from  the  two 
bounded  variation  terms.  This  leads  to  a  proportionality  relationship  between  (3i5  and  p2- 

s 

Since  BV((p5)  =  BV((p)  and  \<j>s  -  </>  |2  =  — ,  and  if  we  take  P2  °=  82  then  we  have 


[14] 


for  some  constant  C.  That  is,  the  three  terms  in  the  cost  functional  are  balanced.  This 
provides  us  with  an  idea  of  values  we  can  use  for  the  parameters.  Pi  and  p2;  namely  that 
the  three  terms  in  the  cost  functional  must  be  roughly  proportional  to  each  other. 
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Mathematical  Formulation 


Now  we  turn  again  to  the  variational  problem.  We  continue  in  the  one¬ 
dimensional  case  as  before.  We  state  the  following  theorem. 

Theorem  1 :  There  exists  a  u  that  minimizes  J. 

Proof.  We  define  a  normed  space,  X,  as  the  space  of  all  continuous  functions  on  (0,1) 
with  bounded  variations,  BV(u)  +  BV(ux)  <  X  is  equipped  with  the  norm 

i 

|*4  =  \\u\dx  +  BV(u )  +  BV(ux ) .  [15] 

o 

Then  it  can  be  shown  that  X  is  complete,  that  is  X  is  a  Banach  Space.  Suppose  {un}  is  a 
minimizing  sequence  of  J.  That  is,  such  thatlim7(Mn)  =  inf  J(u )  over  L^Chl).  Now 

~  n—>* © 

there  exists  a  bound  for  the  subsequence.  Since  X  is  compactly  embedded  in  L  (0,1), 
there  must  exist  a  subsequence  {  u. }  of  {un}  that  converges  strongly  in  L2(0,1)  and 
weakly  star  in  X  to  u.  Since  the  norm  is  weakly  star,  lower  semi-continuous,  we  have 

J(u)  <  liminf  J(u-).  [16] 

n— >oo 

Therefore,  we  have  a  u  that  minimizes  the  cost  functional,  J. 

We  turn  the  discussion  to  the  necessary  and  sufficient  condition  for  optimality. 

To  do  this  first  let  u  be  the  minimizer  for  the  cost  functional  J.  Then  for  all  0  <  t  <  1,  we 
have  J(u  +  t(v  -  u ))  >  J(u )  for  all  v.  Now  we  have 


1  t 2  1 

J  (u  +  t(v  -  u))  -  J  («)  =  j  t(u- z)(v -u)- I — ( v-u)2dx  +  /?jJ|  ux  +t{v -u)x\-\ux\ix 


i 

+  02j\uxx+t(y-u)xx\-\uxx}fix 


Since  we  know  that  for  A,  B  e  R, 
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U  +  t(B  -  A)\  =  1(1  -  t)A  +  tB\  <  (1  -  ON  +  t\B\ 

|  A + »(B  -  A)|  -  |A|  <  i(|b|  -  |a|)  1181 

Using  this  fact,  coupled  with  the  expressing  for  J  above,  we  have 

1  f2  1 

0  <  J(u  +  t(v  -  u ))  -  J {u)  <  J/(w  -  z)(v  -u)  +  —  (v  -  u)2dx  +  tPx  Jjv^l  -  \ux\dx 

0  ,  2  °  [W] 

0 

Now  we  have  the  one-dimensional  case  for  optimality  by  letting  t  — >  0: 

(«-z,v-«)  +  A(|vJt|-|MJ[|)+/?2(|v„|-|«„|)^0  for  all  v  6  L^O.l).  [20] 

Now  let  W={//2(0,l)n^(0,l);  ue  L2( 0,1);  KI)BaeL2( 0,1);  u(0)  =  u(l)  =  o}. 

Then  if  we  assume  that  the  minimizer,  u(x),  is  an  element  of  W.  Then  we  can  derive  the 
differential  form  of  the  necessary  and  sufficient  condition  for  optimality.  First  we  note 
that  J(u  +  th)  >  J(u )  for  all  h  6  W  and  t  e  R.  In  order  to  write  the  optimality  condition 
in  an  equation  form,  we  consider  for  8  >  0, 

7£(M)=ijiM-zi2^+Aiv^-Kr + +K*r  •  [2ij 
^  0  0  0 


As  we  will  see  in  later  discussion,  adding  the  e  avoids  the  difficulty  of  dealing  with 
singularities  like  p-r ,  where  ux  =  0.  Now  doing  the  following  calculation 


■Js2  +| («  +  t/z)J2  -t-^l2 

K«+^r-kr 

-^£2  +  |(m  +  r/z)c|  +  2  + 

2faA+r2|/*,|2 

■^£ 2  +  |(m  +  t/l  )x  +  -^£2  +  Ux 


[22] 
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and  similarly 


iJ72  +|(m  +  ^)J  ~^£2+\uxx\ 

l(«+^LI2-kl2 

iJe2  +  |(w  +  r/i)J2  +^/e2  +|wxr 

_ _ 2fMja./rcc+i  1^1 _ 

+V6-2 +k* 


Applying  these  two  formulas  to  the  optimality  condition  and  forming  a  difference 
quotient,  we  are  left  with 


J(u  +  th)-J(u) 


2 uxhx  + 1 ) 


>J(K-Z)fc  +  -|fc|2G&C  +  AJ-F^- - —  1 

0  +  |(w  +  tf*)J  ~i~yS  +  Ux 


=dx 


[24] 


+AJ 


2m  /i  |2 

a  xr  |  xr  | 


+Vf2 +k 


-dx 


for  all  h  g  W  and  t  e  R. 


Now  taking  the  limit  of  this  expression  as  t  approaches  0,  we  have 

](u-z)h  +  0 I  2+/?2  jMs  =dx  =  0,  [25] 

o  ^e2 +\ux  |  -\/f2+|Mx*| 

for  all  h  G  W.  This  is  the  weak  variational  form.  This  optimality  condition  can  also  be 
written  in  the  differential  form  as 


(u-z)~P\  — 

OX 


V?+K 


+  A 


2&t2 


+K 


=  0, 


[26] 


with  u(0)  =  u(l)  and  uxx(0)  =  uxx(l)  =  0,  for  e  >  0. 


11 


Numerical  Algorithm 

Now  that  we  have  the  mathematics  behind  the  approach  to  the  solution  we  must 
develop  a  numerical  algorithm  that  can  be  readily  implemented.  Using  the  formula  we 
derived  for  the  one-dimensional  case,  we  can  obtain  the  two-dimensional  form  in 
£2  =  (0,1)  x  (0,1)  and  add  a  strategically  placed  8  >  0.  Then  the  formula 

f  ^  f  ^ 


(u-z,<f>)  +  px 


Vw 


■\j  £2  +  Vw 


+  A 


A  u 


tJs2  +  |Aw|‘ 


=  0 


[27] 


is  in  the  differential  form  where  Vu  =  (ux,u  ) ,  the  gradient,  and  Au  =  uxx  +  u  ,  the 


Laplacian.  In  fact,  if  we  define  (p(s)  =  ^£2  +  s  ,  where  s  >  0,  then  it  can  be  shown  that  it 


is  the  optimality  condition  for  the  minimization  of 

Je(u )  =  —  J  |w  —  z|2  +  ^)|vm|2)+  ^(jAM|2)|7x . 


[28] 


We  propose  that  the  fixed  point  iterate  is 


uk+]  -pM 


Vw* 


i 


e2  +  Vw* 


+  P2a 


Aw* 


I 


'•A 


e  +  Aw 


+  juAA(um  -uk)  =  0,  [29] 


where  p,  is  arbitrarily  small  and  is  used  to  insure  that  uk  is  in  the  H2(£2). 

We  state  the  following  theorem  of  the  convergence  properties  of  the  proposed 
algorithm. 

Theorem  2:  The  cost,  Je(uk),  is  monotonically  decreasing  as  the  iteration  progresses.  That 
is  to  say  that  JE(uk+1)  <  Je(uk)  for  all  k.  The  next  part  of  the  theorem  states  that 


lim  J  (u  )  =  inf  J  (u) .  Finally  suppose  that  uk  is  bounded  in  H2(£2).  Then 

u6H2(0,l)nf/J(0,l) 


12 


uk  converges  weakly  to  the  unique  solution  u  in  H2(£2).  That  is  to  say  that  for  all 
(f>  G  H 2  (0,1)  n  Hl0  (0,1) ,  we  have  the  equation 

#(p/(vh|2)vm,V(z>)+  A  (^'(|Vm|2  )^m,  A<z>)+  {u-z,<p)  =  0 .  [30] 
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Implementation  of  Algorithm 

In  our  implementation,  we  will  take  the  proposed  algorithm  and  create  the 
discretized  formulation.  First,  we  establish  that  the  domain  is  the  unit  square, 

(0,1)  x  (0,1)  in  R2.  We  create  a  uniform  n  by  n  size  mesh  of  the  domain  thereby 
2  1 

indicating  n  pixels.  Now  letting  h=  —  ,  and  creating  index  variables  i  and  j  such  that  0  < 

n 

i,  j  <  1.  Then  the  solution  for  u  at  (i,  j),  as  approximated  by  uy,  is  at  the  mid-point  Xi,  j, 


' 

v 

n 

(  l) 

\ 

where  xtJ  = 

i  - 

2  , 

h, 

h 

A 

) 

\  J 

j 

Now  we  represent  the  Laplacian,  -A,  in  terms  of  the  central  difference  scheme  and 
we  obtain  the  matrix 


Sh  =I®H0  +  I®H0,  [31] 

where  ®  is  the  Kronecker  product  and  Ho  is  the  n  by  n  tri-diagonal  matrix  represented  as 


#o  =*‘ 


1  -1 

-1  2  -1 

-1  \ 

2  -1 

-1  1 


Now  we  can  approximate  the  value  of  (p(  Vu  |  2)  at  Xjj  by 


<P\ 


fl[ 

UMJ  -  UiJ 

2 

+ 

uu  ~ui-u 

2\ 

1 

+  — 
2 

V 

2 

+ 

2  V 

2 

v 

h 

h 

J 

h 

h 

)) 

where  the  following  relation  holds 


Un+lJ-U«J  _  Ui»+l  _  Ul<j  ~  U0<j  _  “M  ~“t.O  =q 


h  h  h 

Let  us  now  define  the  following  difference  operators: 


h 


[32] 


,  [33] 


[34] 
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[35] 


Dj+  =  D+  ®  I 

d;  =  d~®i 

D2+  -  /  ®  D+ 
D~  =  l®  D~ 


which  are  in  order  the  forward  difference  in  the  first  variable,  the  backward  difference  in 
the  first  variable,  the  forward  difference  in  the  second  variable,  and  the  backward 
difference  in  the  second  variable.  Now  we  can  give  a  discrete  version  of  the 
regularization  problem  given  by 


AK)=h  -z*f +EE 

;=i  j=i 


|(A  \Uh  +  ^P\  )jM/l|  + |(^2  )jUh  +  iPl  )jUh 


■  fl2tp((Shuh  )!  J )  ] 


The  necessary  and  sufficient  optimality  condition  for  J(u)  is 
ut  +  p\(d;)  a,d;  +(d;) \,d; +(d;) \,Dt  +(d!-)'a,z>2-)1,  +  fi2S„A2Stut  =  z„,  [37] 


2  0  0 
where  zh  e  R"  is  defined  by  (zh)iHj_l)rl  =  z(xI>; )  and  Ai  and  A2  represent  the  n  by  n 


diagonal  matrices  with  the  diagonals  defined  by 


and  Mi+U-i)n  =  9%Shuh)itJ)  respectively. 
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Matlab  Implementation 


To  implement,  we  initially  use  Matlab  to  solve  small-scale  problems  in  their 
entirety.  For  a  small  system,  n  =  70,  we  can  first  create  an  image  (70  by  70  pixels).  Then 
by  adding  a  random  variable  to  this  image,  create  a  noise  element  to  produce  a  noisy 
initial  image.  Then  using  the  algorithm  developed,  we  can  obtain  a  solution.  In  this  case, 
it  is  feasible  to  use  Matlab ’s  built  in  commands  to  solve  the  linear  system  since  it  will  do 
it  relatively  quickly.  Note  that  a  70  by  70  image  results  in  the  need  to  solve  a  linear 
system  of  4761  equations.  That  is  in  general,  an  n  x  m  image  will  result  in  a  linear 
system  of  ( n  - 1  )(m  - 1)  equations.  However  a  70  by  70  pixels  is  a  relatively  small-scale 
image,  and  our  goal  is  to  develop  an  algorithm  that  can  solve  a  256  by  256  problem  in 
minimum  time. 

The  following  is  the  code  used  for  the  main  Matlab  calling  routine. 

% function  [x, e] =reconst (idx, n, H,mu, g,  c , ep, del , bt ) 

%  The  first  part  is  an  initialization  phase  (idx==0) 

%  This  first  block  creates  an  initial  image 
if  idx==0;  m=n+l;  l=n-l;  n0=l*l;  nn=2*n0; 
dx=l/50 ;  tmp= [dx: dx: 1-dx] ' *ones (1,49) ; 

z=min (abs ( tmp) , abs ( 1-tmp) ) ;  z=min ( z , min (abs ( tmp ' ) , abs ( 1-tmp ' ) ) ) ; 
z=min ( z , . 3 ) ; 

zz  =  z;  z  =  zeros (1,1);  z  (10 : 58 , 10:58) =zz;  z ( : , 1 : 25) =zeros (1,25)  ; 

%  The  next  step  adds  in  random  noise  to  the  image  resulting  in  a  noisy 
image 

f0=z+2*del* (rand(l, 1) - . 5*ones (1,1));  f 0=f 0 ( : ) ;  u=zeros (m,m) ; 

%  The  next  step  is  to  create  the  2-D  laplacian 
d=ones (1,1) ; 

h0=n*n*spdiags ( [ -d  2*d  -d] , -1 : 1, 1, 1) ; 

h0=kron ( speye ( 1 ) , hO) +kron (hO , speye ( 1) ) ;  bl=n*spdiags ( [ -d  d] , [ -1 

0 ] ,1,1) ; 

%  Here  we  create  the  forward  and  backward  differences  with  respect  to 
each  coordinate 

b2=kron (bl, speye (1) ) ;  %  Backward  WRT  2nd  coordinate 

bl=kron ( speye ( 1) ,bl) ;  %  Backward  WRT  1st  coordinate 

cl—  bl ' ;  %  Forward  WRT  1st  coordinate 
c2=-b2 ' ;  %  Forward  WRT  2nd  coordinate 
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bb= [bl ; b2 ] ;  cc= [cl ; c2 ] ; 

h^speye (nO ) +mu*hO ;  %  This  line  can  be  omitted  when  solving  large 
systems 

%  x=h\ f 0 ; 

x=  f 0 ;  %  The  above  line  is  replace  with  this  line  to  bypass  the  Matlab 

matrix 

%  solve.  For  larger  systems,  the  matrix  solve  step  is 
computationally 

%  intense.  This  is  valid  since  h  is  approx,  the  identity  matrix 
end  %  End  of  Initialization  Part 
%  Begin  solution  phase  (idx==l) 

%  The  next  step  is  to  calculate  the  sqrt (abs (grad  u)A2) 
if  idx>0; 

ql=bl*x;  q2=cl*x;  q3=b2*x;  q4=c2*x; 

q= . 5* (ql . A2+q2 . A2+q3 . A2+q4 .  A2 ) ;  kl=f ind (q< . 2 ) ;  k3=f ind (q>l ) ; 

%  This  next  s 

if  idx==l;  qq=sqrt (ep+q) ;  b=g./qq;  end 

a=spdiags ( [b; b]  ,  0 , nn,nn) ; 

c=gg/n/n;  qq=hO*x;  qq=c . /sqrt (ep+qq. A2 ) ; 
aa=spdiags (qq, 0 , nO , nO) ; 

%  The  next  step  is  to  find  the  hessian  matrix 
hes=hO*aa*hO+bb' *a*bb+cc ' *a*cc+h;  f =f 0~hes*x; 

%  The  following  is  the  interface  for  calling  the  multigrid  solver 
n 

(n-1) A2 
nnz (hes) 

dummy= input (' Pausing  ...  edit  for  new  n?  1  yes,  2  no:  '); 
if  dummy— 1 

unix('emacs  test.f  &'); 

end 

dummy= input ( x Pausing  ...  recompile  for  new  n'); 
if  dummy==l 

unix('f77  test.f  amglr5.o  ctime.o7); 
end 

%flag=l 

test0007 

%flag=2 

unix ( 7 a. out ' )  %  This  is  the  actual  solve  step. . .replaces  matlab  matrix 
solve 

dx=load( 7ansu7 ) ;  x=x+dx; 

u(2 :n, 2 :n) =reshape (x, 1, 1) ;  end  %  Solution  Phase  (idx==l) 
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Multigrid 

Our  problem  of  speed  reduces  to  the  age-old  problem  of  solving  a  linear  system. 
In  addition  to  getting  a  robust  and  accurate  answer,  speed  is  now  an  important  factor.  We 
turn  to  the  use  an  algebraic  multigrid  method  to  solving  the  system.  Multigrid  is  a  coarse 
grid,  relaxation  method.  In  multigrid  operations,  there  are  typically  four  steps  to  the 
algorithm.  The  first  involves  building  a  coarse  grid  of  the  original  system  and 
performing  relaxation  on  a  fine  grid  until  the  error  becomes  smooth.  The  relaxation  step 
tends  to  smooth  errors.  This  dampening  effect  tends  towards  more  accurate  solutions. 
After  computing  and  then  restricting  the  residual,  it  is  transferred  to  the  coarse  grid. 
Solving  the  coarse  grid  residual  equation  is  the  next  step.  Then  there  is  the  interpolation 
of  the  error  and  the  correction  stage. 

Algebraic  multigrid  is  a  more  robust  algorithm  that  does  not  take  advantage  of  the 
geometry  of  a  system.  Here  one  can  apply  multigrid  to  an  unstructured  grid,  which  is 
impossible  to  do  in  the  geometric  multigrid  scheme.  In  essence,  it  is  the  opposite  of 
geometric  multigrid.  First  the  algebraic  method  defines  the  multigrid  components  and 
then  it  performs  the  multigrid  cycles.  More  specifically  it  fixes  the  relaxation  and  then 
defines  the  multigrid  components  so  that  coarse  grid  correction  eliminates  error  not 
reduced  by  relaxations. 

There  are  two  phases  to  algebraic  multigrid.  The  first  is  the  setup  phase.  It 
involves  setting  up  the  coarse  grids  and  defining  interpolation,  restriction,  and  the  coarse 
grid  operators.  The  second  step  is  the  solution  phase.  Here  we  have  the  standard 
multigrid  cycling  operations. 
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Fortran  AMG 


Now  that  we  have  described  a  little  about  how  algebraic  multigrid  works,  we  must 
now  put  it  into  practice.  First  I  obtained  a  recommended  algebraic  multigrid  package 
from  http://www.mgnet.org/mgnet-codes-gmd.html.  This  package  is  a  Fortran-77  driven 
multigrid  package  developed  by  John  Ruge,  Klaus  Stueben,  and  Rolf  Hempel. 

The  package  we  used  contained  four  Fortran  files:  amglr5.f,  auxlr5.f,  ctime.f, 
drvlr5.f.  The  first  file  contained  code  that  was  the  entire  algebraic  multigrid  process. 

The  second  file  was  a  driving  interface  program  that  setup  parameters,  input,  and  output 
methods.  The  last  two  files  are  not  used  since  they  do  not  pertain  directly  to  the  solution 
of  the  system.  However  ctime.f  must  be  compiled  and  included  when  creating  the  main 
executable  file.  drvlr5.f  does  not  need  to  be  included.  In  using  this  package,  we  did  not 
use  several  available  features  for  ease  of  use.  The  main  subroutine,  amglr5.f  was  not 
changed.  The  driving  program  auxlr5.f  was  changed  to  fit  our  needs.  Before  we  get  into 
change,  we  must  first  describe  how  information  is  stored  so  that  it  is  usable  to  this 
program. 

Given  a  linear  system  L  ■  x  =  /  ,  we  need  to  store  this  information  in  a  series  of 
vectors.  First  there  are  several  assumptions  on  L.  The  first  is  that  diagonal  elements  are 
always  positive  on  all  grids.  The  second  requirement  of  the  program  is  that  L  is  a  square 
matrix,  which  is  singular  with  row  sums  equal  to  zero.  These  are  requirements  of  the 
program.  Additionally  there  are  some  theoretical  restrictions  to  L.  L  is  a  positive 
definite,  most  of  the  off  diagonal  elements  are  nonpositive,  and  row  sums  are  non¬ 
negative.  Since  we  are  working  with  large,  sparse  matrices,  we  need  to  develop  an 
efficient  storage  scheme.  The  information  given  by  the  matrix  is  stored  in  a  series  of 
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three  vectors  in  what  is  called  “compressed  sky-line”  fashion.  The  first,  calling  it  A, 
contains  the  non  zero  elements  of  L  stored  by  rows,  with  each  row  starting  with  its 
diagonal  element,  and  the  rest  of  the  other  non-zero  entries  following  the  leading 
diagonal  entry  in  any  order.  The  last  entry  of  A  contains  a  dummy  element,  which  is 
needed  to  match  the  sizes  of  the  other  vectors.  It  contains  no  usable  information  and  its 
value  is  arbitrary. 

Now  in  order  to  preserve  the  locations  of  the  non-zero  elements  in  the  sparse 
array,  we  create  two  additional  vectors,  IA  and  JA.  Now  letting  NNU  represent  the 
number  of  unknowns  in  the  system  and  NNA  represent  the  total  number  of  non-zero 
elements  in  A,  we  can  define  IA.  Elements  of  IA(I)  point  to  the  position  of  the  diagonal 
entry  of  row  I  with  in  the  vector  A.  For  example  if  the  third  diagonal  element  is  stored  in 
the  6th  entry  in  vector  A,  then  IA(3)  contains  the  value  6.  Then  theoretically  the 
dimensional  of  vector  A  should  be  NNU.  However  the  program  needs  the  length  of 
vector  A  to  be  NNU+1,  with  the  last  entry  being  NNA+1.  This  last  addition  allows  the 
program  to  recognize  the  length  of  the  last  matrix  row.  To  define  JA,  we  first  establish 
that  is  has  the  same  length  as  A.  That  is  it  contains  NNA+1  elements.  Remember,  when 
we  defined  A,  there  was  an  extra  dummy  element  attached  to  the  end.  Now  the  Jth 
element  in  JA  corresponds  to  the  Jth  element  in  A.  More  specifically,  JA(J)  points  to  the 
column  the  entry  A(J)  was  in,  in  the  original  matrix  L.  For  example  if  the  (4,5)  element 
of  matrix  L  was  -3,  and  that  entry  was  the  9th  entry  in  vector  A,  then  JA(9)=5.  Here 
again,  since  the  length  of  JA  must  match  that  of  A,  then  another  arbitrary  value  must  be 
added  for  the  NNU+1  entry  of  JA. 
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Now  as  it  seems,  the  storage  requirements  of  these  vectors  are  known  ahead  of 
time.  However,  since  the  routine  will  actually  change  these  vectors  in  the  multigrid 
process,  we  need  to  allocate  much  more  space  than  these  vectors  initially  use.  The 
program  suggests  that  the  size  of  the  vectors  A  and  JA  should  be  3*NNA  +  5  *  NNU 
while  the  size  of  IA,  U,  and  F  should  be  2.2  *  NNU.  And  the  size  of  IG  is  5.4  *  NNU.  In 
practice  however,  we  found  that  IG  would  often  need  more  allocated  memory,  so  10  * 
NNU  was  used  instead.  Vector  F  contains  the  “right  hand  side”  of  the  system.  Output 
vectors  U,  which  contains  the  solution,  and  IG,  which  is  used  internally,  are  two  vectors 
which  space  must  be  allocated  to  beforehand. 

Now  we  can  look  at  the  call  to  the  main  subroutine.  In  the  code,  the  line  is 

CALL  AMG1R5 (A, IA, JA, U, F , IG, NDA, NDIA, NDJA, NDU, NDF, NDIG, 

+  NNU, MATRIX, ISWTCH, IOUT, I PRINT, LEVELX, 

+  IFIRST, NCYC , EPS , MADAPT, NRD, NSOLCO, NRU, 

+  ECG1 , ECG2 , EWT2 , NWT, NTR, IERR) 

The  first  6  parameters  have  been  mentioned  previously.  The  next  six  parameters  NDA, 
NDIA,  NDJA,  NDU,  NDF,  NDIG  represent  the  dimensioning  of  the  corresponding 
vector  in  the  calling  program.  NNU  is  again  the  number  of  unknowns.  MATRIX  is  a 
two-digit  integer  variable  that  represents  what  kind  of  matrix  L  is.  The  first  digit  is  a  1  if 
L  is  symmetric  and  2  if  it  is  not  symmetric.  The  second  digit  is  1  if  L  has  row  sum  zero 
and  2  if  it  is  not.  In  our  use,  MATRIX  will  have  a  value  of  12,  a  symmetric  matrix 
without  rowsum  equal  to  zero. 

The  rest  of  the  parameters  are  used  to  control  how  multigrid  works.  EPS  is  the 
value  for  convergence  criterion.  This  value  was  most  critical  in  finding  out  the 
correlation  between  accuracy  and  calculation  times  of  the  program.  ISWITCH  is  a 
parameter  controlling  which  modules  of  AMG1R5  are  being  used.  Here  we  use  a  value 
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of  4,  which  allows  for  all  four  modules  to  be  used.  The  first  module  is  SETUP,  which 


defines  the  operators  needed  in  the  solution  phase.  The  second  is  FIRST  which  initializes 
the  solution  vector.  SOLVE  computes  the  solution  by  algebraic  multigrid  cycling.  And 
the  final  module,  WRKCNT,  provides  information  about  residuals,  storage  requirements, 
and  cycles. 

IFIRST  is  a  value  that  represents  the  parameter  for  the  first  approximation.  It  is 
also  a  two-digit  number.  The  first  digit  is  arbitrary  since  it  is  not  used  and  is  set  to  1, 
since  it  must  be  non-zero.  The  second  digit  is  0,  1,  2,  or  3.  0  represents  no  setting  of 
first  approximation.  1  represents  the  first  approximations  of  zero.  1  represents  a  first 
approximation  of  1.  3  represents  the  first  approximation  to  be  a  random  function 
determined  by  the  program. 

IOUT  is  a  parameter  that  represents  what  output  is  shown  to  the  user.  For  our 
purposes  we  set  this  integer  value  to  be  13,  which  specifies  residual  output  information. 
BPRINT  is  a  five-digit  integer  value  that  represents  Fortran  unit  parameters.  We  set  this  to 
be  10606.  IERR  is  an  output  parameter,  which  contains  error  flags.  A  value  of  zero 
translates  to  no  errors.  A  negative  value  will  be  a  non-fatal  error,  while  a  positive  value 
represents  fatal  errors. 

LEVELX,  NCYC,  MADAPT,NRD,  ,NRU,  ECG1,  ECG2,  EWT2,  NWT,  NTR 
are  parameters  which  have  standard  values.  LEVELX,  MADAPT,NRD,  ,NRU,  NTR 
have  standard  values  of  zero.  NCYC  has  a  value  of  10250.  NWT  is  2.  ECG1  is  0.D0, 
ECG2  is  0.25D0,  and  EWT2=0.35D0.  Finally  NSOLCO  is  a  parameter  that  toggles  the 
use  of  an  external  package  not  provided  by  AMG1R5.  We  set  this  parameter  equal  to  1 
to  turn  this  feature  off. 
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Results 


In  this  section,  we  describe  how  the  proposed  method  performs  using  a  simulated 
test  example.  Many  factors  can  come  into  testing  our  proposed  algorithm.  The  first  is 
what  are  the  performance  concerns.  Obviously  there  is  a  degree  of  accuracy  involved.  In 
these  types  of  inverse  problems,  accuracy  is  often  hard  to  determine.  But  since  our 
testing  involves  an  initial  image  that  is  polluted  with  a  certain  amount  of  noise,  we  can 
have  a  measure  of  accuracy.  First  we  create  the  test  image,  zy,  as  shown  in  Fig.  4. 
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Fig.  4 

Then  we  add  the  noise  to  the  test  image  resulting  in  the  observed  image,  yy. 

yi,j=ziJ+SniJ,  [39] 

where  ny  is  an  independent,  identically  distributed  random  variable  with  uniform 
distribution  on  (-1,1), .  This  noisy  image  is  shown  in  Fig.  5. 
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Fig.  5 

Accuracy  in  large  part  is  going  to  be  determined  by  our  choices  of  constants,  Pi 
and  ^2-  These  two,  as  mentioned  before  in  the  mathematical  discussion,  control  the  edge 
preserving  nature  of  the  method  and  the  shape  preserving  method.  Also  a  concern  is  the 
phenomena  of  the  deterioration  of  the  height  (level  of  gray  scale)  of  the  image.  Finding 
an  acceptable  balance  that  keeps  edges,  smoothes  slopes,  and  does  not  decrease  the 
height  of  the  image  is  our  goal. 

In  doing  so  it  is  necessary  to  test  different  values  of  Pi  and  P2  together.  What  we 
observed  is  that  a  smaller  value  of  Pi  will  result  in  less  of  the  deterioration  phenomena. 
This  effect  is  shown  in  Figs.  6,  7,  and  8.  For  each  of  these  examples,  P2  =  0.  The  values 
for  Pi  are  0.01, 0.02,  and  0.1,  respectively.  Notice  that  the  larger  Pi  is,  the  blockiness  is 
accentuated.  Also  as  Pi  increases  in  value,  the  deterioration  effect  increases. 
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We  now  introduce  the  second  BY  term.  Here  we  are  careful  to  pick  our  value  of  P2  so 


that  it  removes  the  blockiness  in  the  smooth  regions  and  retards  the  deterioration  effect. 
Fig.  9  represents  the  BV  regulation  method  with  choice  of  P2  =0.007.  The  values  for  P2 
are  0.01,  0.05,  0.1  used  successively. 


Fig.  9 

The  parameter  values  for  Fig.  10  are  Pi  =0.01.  The  values  for  P2  are  0.01,  0.05,  0.1  used 


successively. 
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As  one  can  see,  careful  choice  of  parameters  will  result  in  a  reconstructed  image  with  the 
three  types  of  image  composition.  Here  we  see  that  increasing  Pi  increased  edge 
definition,  but  the  successive  choices  for  p2  allowed  us  to  retard  blockiness  and  the  height 
deterioration  effect. 

Another  concern  is  the  speed  in  which  an  answer  is  produced.  Because  the  size  of 
the  system  increases  with  the  square  of  the  number  of  pixels,  large  images  can  result  in 
time-consuming  calculation.  Our  method  using  multigrid  was  chosen  to  decrease 
calculation  times,  however  careful  choice  of  epsilon  may  also  create  significant  gains  in 
compute  times.  An  epsilon  10  6  vs.  10'14  produced  the  same  accuracy.  The  number  of 
AMG  cycles  performed  was  less  when  we  used.  This  performance  is  based  on  a  128  by 
128  pixel  image.  Testing  of  256  by  256  was  interrupted  due  to  storage  limitations.  A 
better  interface  between  the  Fortran  AMG  routines  and  the  Matlab  code  would  remedy 
this  problem. 
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